#include <stdio.h>
#include "mpi.h"
#include <math.h>
#include <stdlib.h>

/* Global values and subroutines */
#define MAXN 1000000
double x[MAXN];   // Positions of particles
double move();    // Routine for changing position of particles
                  // based on current position and applied force

void update();
void updateNN();
double averageDistance();
void init();


/* Other variables */
int N = 1000;         // Number of particles (N<MAXN)
int Ntime = 100;      // Number of timesteps
double deltaT = 0.1;  // Timestep

int numprocs = 0;     // the total number of processors
int procnum = 0;      
int root = 0;         // the root processor

int main(int argc,char* argv[])
{

   int i, randseed;
   double avgeDist;

    MPI_Status status; /* Return status for receive */

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &procnum);
    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);


   if(procnum == 0)
   {
   	   /* Read in parameters */
      printf("Enter number of particles [e.g. 1000]: ");
      scanf("%d", &N);
      printf("Enter number of timesteps [e.g. 100]: ");
      scanf("%d", &Ntime);
    }

   // Broadcast the timesteps and num of particles
   MPI_Bcast(&Ntime, 1, MPI_INT, 0, MPI_COMM_WORLD);
   MPI_Bcast(&N, 1, MPI_INT, root, MPI_COMM_WORLD);

   /* Initialize positions */
   init();

   /* Print out average position */
   avgeDist = averageDistance();
   if(procnum == 0)
   printf("0 Average distance between particles at timestep 0 is %f\n\n",
          avgeDist);

   /* Loop over timesteps for update using all particles */
   if(procnum == 0)
   	   printf("Using update based on forces over all particles\n");
   for (i=0;i<Ntime;i++)
   {
      update();
   }

   /* Print out average position */
    avgeDist = averageDistance();
    if(procnum == 0)
        printf("update Average distance between particles at timestep %d is %f\n\n",
           Ntime, avgeDist);

   /* Reinitialize positions to the same values as before */
   init();

   /* Loop over timesteps for update using nearest neighbour particles */
   if(procnum == 0)
       printf("Using update based on forces from nearest neighbour particles\n");
   for (i=0;i<Ntime;i++)
   {
     updateNN();
   }

   /* Print out average position */
   avgeDist = averageDistance();
   if(procnum == 0)
   printf("updateNN distance between particles at timestep %d is %f\n\n",
          Ntime, avgeDist);
   
   // End of the program 
   MPI_Finalize();
   return 0;
}



/**************************************************************/
/* Initialize the position values to be evenly distributed.   */
/* (alternatively could give them random initial values)      */
/**************************************************************/
void init() {
   int i;
   for (i=0;i<N;i++) {
      x[i] = (double)i / (double) N;
   }
}


/**************************************************************/
/* Update position of particle using a bogus version of       */
/* Newton's laws of motion (ignore masses and velocities).    */
/**************************************************************/
double move(double xinit, double force)
{
   double xfinal;
   xfinal = xinit + 0.5 * force * deltaT * deltaT;
   return xfinal;
}
